Elastic Instability Triggered Pattern Formation 
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Recent experiments have exploited elastic instabilities in membranes to create complex patterns. 
However, the rational design of such structures poses many challenges, as they are products of 
nonlinear elastic behavior. We pose a simple model for determining the orientational order of 
such patterns using only linear elasticity theory which correctly predicts the outcomes of several 
experiments. Each element of the pattern is modeled by a "dislocation dipole" located at a point 
on a lattice, which then interacts elastically with all other dipoles in the system. We explicitly 
consider a membrane with a square lattice of circular holes under uniform compression and examine 
the changes in morphology as it is allowed to relax in a specified direction. 



I. INTRODUCTION 

It is a testament to the ingenuity of Nature that at 
all length scales there exists a multitude of complex self- 
assembled patterns. The formation of patterns, from the 
textured dimples in a grain of pollen to the labyrinth of 
creases in a coral colony to the ridges in a fingerprint, 
is consistent and repeatable [1-3]; the mechanisms that 
drive them are far from understood. Mastery of such 
processes would revolutionize the fabrication and design 
of novel materials with specific properties. By harness- 
ing elastic instabilities in elastomeric membranes [4] , the 
once lofty goal of creating self-assembled complex pat- 
terns with long range order is now a one step closer. 

When an elastomeric membrane with a square lattice 
of circular holes (with diameter roughly half the lattice 
spacing) in it is uniformly swollen or, alternatively, is 
compressed hydrostatically, the holes deform into a dia- 
mond plate pattern. The order persists for upwards of 
10^ times the original lattice spacing with only phase-slip 
defects, which do not change the overall symmetry of the 
pattern. It is interesting to note that the same diamond 
plate pattern emerges from membranes with vastly differ- 
ent lattice spacings, ranging from 1 fim to 1 cm. Because 
the same mechanism causes patterns to form over such a 
wide range of length scales, this technique is well suited 
to create devices for many fields and industries. In par- 
ticular, when a material is in the diamond plate state, it 
has a photonic band-gap [5]; hence, merely by compress- 
ing and relaxing a membrane, we can control its band 
structure. 

This effect occurs in the highly non-linear regime of 
elasticity; thus, only finite element simulations using spe- 
cific models of nonlinear elasticity capture the entire pro- 
cess of the holes collapsing [6]. However, some simula- 
tions merely predict the orientational order of the col- 
lapsed holes and cannot capture the details of their final 
shape. With the sole assumption that each hole collapses 
to some elongated shape, our model uses only linear elas- 
ticity to successfully predict the orientational order in 
the diamond plate pattern and the herringbone pattern 
formed from an underlying triangular lattice. Not only 
does our model shed light on the interactions in the sys- 



tem, it greatly facilitates the rational design of other pat- 
terns and devices. 

In the next section, we review the linear theory of elas- 
ticity and consider a hydrostatically compressed mem- 
brane containing either a single hole or two holes. The 
breakdown of linear elasticity illustrates the need for non- 
linear analysis. In section III, we adopt the spirit of the 
theory of cracks and show that complex behavior can 
be explained by recasting the problem using linear dislo- 
cation theory. While an explicit analytic description of 
the shape of a collapsed hole cannot be reached through 
this method, each hole can be modeled as an elastically 
interacting distribution of parallel edge dislocations, or 
"dislocation dipoles," whose centers are fixed to the cen- 
ter of their corresponding hole but are allowed to rotate 
freely. In the simple case of a quartet of holes, the result- 
ing configuration exactly reproduces the unit cell for the 
diamond plate pattern. Since elastic interactions are long 
ranged, the collective interactions of all the holes must be 
included to find the true groundstate of the system- the 
inclusion of which additionally stabilize the pattern. We 
end with technical appendices which argue the validity 
of our approximation. 



II. IN WHICH THE LINEAR THEORY OF 
ELASTICITY FAILS 

A. Linear Elasticity in a Nutshell 

The theory of linear elasticity describes the deforma- 
tion and energetics of a solid body under external force 
or load. While this theory was originally developed over 
a century ago and has been presented time and again 
[7-10], this brief tutorial will illuminate a few salient fea- 
tures as well as familiarize the reader to our notation. 
The goal is to develop the framework for a generalization 
of Hooke's Law for solid three dimensional bodies. When 
a solid is deformed, the displacement of every point, x 
is described by the vector it, such that its final positions 
are given by a;' = a; + u{x). Let us first consider two 
points separated by ds — \/ dx ■ dx that are very close 
together. After being deformed, their separation becomes 
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(is', where ds'"^ = (dxi + duiY using the Einstein sum- 
mation convention. By noting dui = {dui/dxh)dxk, we 
may rewrite this as 



ds- = a,2^y(^^p^^p 
^ \dxk oxi dxk oxi 

= ds^ + 2 ufi^dxidxk, 



dxidxk 



(1) 



where is the Lagrangian strain tensor. However, for 
small deformations we need only consider terms linear in 
dui/dxk and the linearized strain tensor is 



1 / duj duk \ 



2 \dxk dxi J 



(2) 



In the following, we will rely upon orthogonal coordi- 
nates, ^fj,, with the diagonal metric ds"^ = h'^da^ + 
h'^df}'^ = [-'--'-]• general, the orthogonal 

coordinates are nonholonomic but the benefit of orthog- 
onality outweighs this complication. 

What is the linearized Lagrangian strain tensor in our 
new coordinates? By definition 



(3) 



and 



Since m is a vector, we may write it in either coordinate 
system: 



(5) 



from which we have Ui = u^^^ ■ Xi. In order to cal- 
culate the direction cosines, we note that if Xi = 
then 



Edxi ^ 



(6) 



We recognize the radicand in the denominator as /i^ and 
we have 



dxj 



dxk , 





dxi 


d 




)- 






\hpdip)_ 



d^^. (7) 



Using orthogonality, J^i ff^ff^ — ^It^ttvi is straight- 
forward to find 
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1 dxj \ dx 



h _L h ^}!2L 



1 dxi \ dx 



d 



dS,v \hp dS,pJ 



(8) 



where we only sum over the repeated indices i and p, 
not n or i> associated with the scale factor h. In the fol- 
lowing we will adopt the Einstein summation convention 
and, only when there is ambiguity, will we specify which 
indices are to be implicitly summed. 

In the spirit of Hooke's Law, there is an energy cost as- 
sociated with displacing every point from its equilibrium 
position. While the general form is quite complicated, to 
first order the energy can be constructed from the linear 
strain tensor. Because the strain tensor is a symmetric 
rank two tensor, the only two possible scalar invariants 
that can be constructed arc {uu)"^ and {uikY . Because 
these terms are invariants of the system, they cannot de- 
pend on coordinate system. From this we deduce the 
form of the energy density 



(9) 



where A and ji are Lame coefficients. Since any deforma- 
tion may be written in terms of uniform or hydrostatic 

compression and pure shear, we rewrite the strain tensor, 
Uik = \kkUii + {uik - \5ikUii). The hydrostatic compres- 
sion is given by the first term, as it involves only the 
trace of Uik- And the second term describes pure shear 
because its trace is zero. In terms of these quantities, the 
energy becomes 



(10) 



where the bulk modulus K = X + ^1-^, and ^ is the shear 
modulus. 

When a body is deformed, the displaced internal el- 
ements experience forces which tend to restore them 
to their equilibrium positions. The volume element 
bounded by the surfaces x = xq, x = xq + Sx, y = yo, 
y = ya + 5y, z = zq, and z = zq + (5z experiences a force 
along any surface with normal n, /j = aikUk, where aik 
is the stress tensor and 



(11) 



In equilibrium, the internal stresses of the system must 
balance the external forces exerted upon it; thus, the 
equilibrium condition for the system is 



dxk 



- /r = 0. 



(12) 



In curvilinear coordinates, this entire discussion can be 
repeated and the stress tensor in these coordinates is 
merely a transformation of the stress tensor in Cartesian 
coordinates: 



1 1 



dxi dxk 



K ^ d^p, 



fife, 



(13) 



where, again, there is no sum over jj, or v. The equilib- 
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rium conditions become 



)_dx^ 



1 



+ 



1 do 



- /r = 0- 



(14) 



where we have followed the same procedure as in the 
derivation of u^^ which requires expanding the tensor in 
both the Xi and frames. 

The relation between stress and strain follows the ar- 
gument in Cartesian coordinaties. If the system is de- 
formed an infinitesimal amount 5ui, the work done by 
the change in internal stresses is the force times the dis- 
placement, W = J dV {d(Jik/dxi)5uk- When integrated 
by parts, the work is 



W = (t) dA cFik rik 6ui - j dVaik 
■ j dVaik duik- 



dSui 
dxk 



ddui 
dxk 



(15) 



Note that the surface integral vanishes because 
CTjfc = at infinity. Thus, dE = aik Suik ■ 
By taking the total differential of equation (10), 
dE = Kuiiduii + 2fj.{uik - 2^ikUu)d{uik - ^^ikUu) = 
{Kuu 5ik + '2^J.{ulk - 2^ikUii)) duik, we can rewrite the 
stress tensor in terms of the strain tensor, aik = KSikUii — 
2^* (uik — 2^ikUu) , or conversely, the strain in terms of 
the stress, Uik = -^Sikcra + ^ {aik - ^^ikcrii) ■ 

Similarly, it is instructive to consider a homogeneous 
deformation, which has constant stain tensor everywhere 
in the volume. For a rf-dimensional solid, uniform pres- 
sure is applied to the faces with normals in the ±z 
directions. In 3-dimensions, for example, we consider 
the simple compression of a rod. This implies that 
azin-i = p, or (Tzz = P- Thus, all off diagonal com- 
ponents of the strain tensor are zero, and the diagonal 
components are uu = p{-^ — j^)/^, for all I ^ z, and 
Uzz = + The relative longitudinal compres- 

sion is given by u^z = p/Y^, where 



2d'^Kn 



2ii + (^2 - d)K 



(16) 



is the d-dimensional Young's modulus. The ratio of 
transverse extension to longitudinal compression is given 
by the Poisson ratio 7 = —uu/uzz = {dK — 2ii)/{{d'^ — 
d)K+2fjb). Conversely, the bulk and shear moduli written 
in terms of the Young's modulus and Poisson ratio are, 
respectively, K = ^Yd/ {l-{d-lYi) and /i = \Y,i/(l+-i). 
The stress tensor is given in terms of the strain tensor by, 



fife 



Yd 
1 + 7 



Uik 



l-(d-l)7 



5ikUii 



(17) 



and the converse by, 
1 



Uik = ^ ((1 + ~ l5ikcru) . 



(18) 



The conventional form of the energy is given in terms of 

Yd and 7, 



E = 



Yd 



2(1 + 7) 



d'^x ( Uik^ + 



7 



l-(d-l)7 



uu' . (19) 



Now that we have derived the relations between the stress 
and strain tensors, the equilibrium condition may be re- 
cast in terms of the displacement vector, 



Yd, 



ay l-(d-3)7 d^H 
2(1 + 7) V 54 1 - (d - 1)7 dxidxi 



/. = 0. (20) 



If the force only acts through the surface, then fi vanishes 
in the bulk and we recover the result that V^V • u = 
from which it follows that V^m = 0. 



B. Elasticity in Flatland[12] 

Our ultimate goal is to study the effects of uniform 
tension on a thin sheet of elastic material with a square 
lattice of circular holes cut in it, as it models an elastic 
sheet that is uniformly swollen. In general, the result- 
ing deformations are constant across the thickness of the 
sheet and may be considered to be purely longitudinal. 
This is because forces act primarily in the plane of the 
film, yielding the boundary condition (Jikiik = 0. Hence- 
forth, we will only consider two dimensional systems with 
planar deformations. 

Let us pause for a minute to derive the equilibrium 
conditions of this system. Since the deformations are 
constant throughout the thickness of the sheet, we may 
assume that Uzz = 0, and thus azz = ^iz = 0. We 
are left with the two equilibrium equations hduik /dxk = 
—pi, where i,k = x, y, and h is the thickness of the film. 
When no external body forces are present the equations 
of equilibrium reduce to daik/dxk = or. 



dx 



dy ' dx 



+ 



dy 



= 0. 



(21) 



Viewing these both as equations of the form V • A = 0, 
we know that Uxi = Ukdk4>y ^iid Ojy = tjkdk4>x and it 
follows that dx(f)y + dy(f>x = Cxy ~ ^xy = (where Cik is 
the totally antisymmetric tensor). Thus (j)m = f-mndnX 
for some scalar x- known as the Airy stress fiuiction. 
We thus have aik = ^im^kndmdnX- Moreover, since 
Uik = {diUk + dkUi)/2, we have eimekndidkUmn = 0. The 
relation between stress and strain, (17) implies that 



, ^ Yd 

dxidxk 1+7 



+ 



7 



l-(rf-l)7 



= 
(22) 
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As a result, we deduce that the stress function satisfies 
the biharmonic equation, V^x = 0. In orthogonal coor- 
dinates, we recast the components of the stress tensor in 
terms of the Airy stress function by taking advantage of 
Eq. (13) and reexpressing dij in terms of ct 



E 



1 



_d_ 



1 dhp dx 



dip J ""h^K d^^di. 



(23) 

where = 1 because {a, (3} is a right-handed orthonor- 
mal basis. Equations (23) are solutions to the equilibrium 
equations, Eq. (14), and the Airy stress function solves 
the biharmonic equation in orthogonal coordinates. 



1. Fixing a Hole (Demo) 

The simplest system to study is an infinite elastic sheet 
with a circular hole, of radius R, cut in it under uniform 
tension PSc. This problem naturally lends itself to polar 
coordinates, for which the equations for the stress func- 
tion become, 



Idx 
r dr 
d_ 
dr 



1 d^x 



Idx 
r 86 



d'^X 



A2 ' 



(24) 



A standard procedure for solving such problems is to 

solve first for the deformation of a continuous sheet un- 
der the proper forces. Secondly, we solve a for a second 
stress function respecting the symmetry broken by the 
force with boundary conditions aik{r = oo) = 0. The 
final stress function is given by the sum of the two stress 
functions, where the matching condition is given by the 
stress free boundary condition at the edge of the hole. 

The components of the stress tensor for a continuous 
elastic sheet under uniform tension Px are a^} = P and 



.(0) - ^(0) 



0, which, by integrating Eqs. (21) yield 



the stress function 



X 



Pr^ (sin^ (/)) /2 = Pr^ ( 1 - cos 20) /4 (25) 



from which it follows that the components of the stress 
tensor arc a'i^J (r) — P(l -|- cosi 



)/2, af],{r) = P(l 



cos20)/2, and cr^^ (r) = P(sin20)/2. Clearly, rotational 
symmetry is broken in the x-direction. 



In order for the second stress function to respect the 
broken symmetry of the system, it must have the form 
x'^' = /(^)+5(^) cos 2(1). Since the stress function satisfies 
the biharmonic equation, we can easily integrate to find 



/(r) = ar^ log r -f -|- c log r 
g{r) = sr'^ +tr'^ +u/r'^ + v 



(26) 



The first boundary conditions (T^^^r = 00) = dictate 
that a = b = s = t = 0, leaving 

(7W(r) = c/r^ - (6u/r^ + Av/r'^) cos 20 

= -c/r2 + 6(w/r-^)cos20 
a-^y(r) = - (6u/r* + 4w/r^)sin20 (27) 

Using the final boundary conditions, a^^J (R) = — <tJ:"'(P) 
the remaining constants arc foimd to be c = — PP^/2, 
u = -PR'^/4:, and v = PR^/2. Note that even though 
there are only two equations for three unknowns, this 
system is not underdetermined because the constant c 
cannot depend on 0. Assembling this, the components 
of the stress tensor are given by 



(Trr{r) 



p 
2" 
p 
2" 



1 + 



P2 

R^ 



1 - 



1 + 



4P2 
3P^ 



3P^ 
cos 20 



cos 20 



P 
2" 



1 + 



2r2 



)sm20. 



(28) 



which may be rewritten, using equation (18), as compo- 
nents of the strain tensor 





1 . 
Y2 






-12 


- JCTrr) 


'^r4> — 


1+7 

'■2 





(29) 



Recall, the strain tensor is the relative displacement 
of every element from its equilibrium position. In polar 
coordinates we have Urr = dur/dr. u,p,p = {du^ / d(f)) / r + 
Ur/ r, and 2ur0 = duc/, /dr+ {dur/dcp) — u^/r, from which 
it follows that the displacement vector is 



Ur = 



W0 



p 
2Y2 



^ (1 + 7)^' A, ^/ R"" 



+ 



4P^ 



cos 20 



,2\2 



P (fi^ + r^)^+7(i;^ - r2) 
2Y2 r3 



■ sin 20. 



(30) 
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However, if we repeat the above process for an infinite 
sheet under hydrostatic compression (or expansion), to 
hnear order the rotational symmetry of this system is 
not broken. Consider an annulus of inner radius Rq and 
outer radius i?i under uniform hydrostatic compression 
with boundary conditions given by (Trr' — —P and cr^'^' = 

a^^J — 0. Clearly, the displacements are purely radial, 
and we need only solve V'm = 0, subject to the bound- 
ary conditions arrir = i?i) = —P and arri^ = ^o) = 0- 
The displacements are given by Ur{r) = a/r + br; from 
which it follows arr = [—(1 — ^)a/r'^ + (1 + 7)^] ■ 

The boundary conditions determine the values of the con- 



a = 



stants a ■ 



(1 + 7) 



and b 



The solution to an annulus under hydrostatic compres- 
sion is 



R1 



Urir) 



p 



7?2 



(l+7)^ + (l-7)r 
r 



Y2RI 

with the components of the stress tensor given by 

Rl 



(31) 



.(r) = 



^) = -P 



Rl 



Rl 



Rl 



Rl 



Rl 



Rl 



Rl 



<7r4,ir) = 0. 



(32) 



In the case of a finite sheet or a pipe under hydrostatic 
compression, there is the well known von Mises buckling 
instability at a critical pressure where the circular hole 
deforms into an ellipse whose major axis is chosen at 
random [13]. This critical pressure scales with the ratio 
of the system size to the hole radius and, thus, diverges 
for large systems. The system we are studying, on the 
other hand, has an underlying lattice which breaks the 
rotational symmetry of each hole. A superposition of 
the above solutions would not account for the interaction 
between holes. 



2. Fixing a Hole (Take 2) 

Understanding the elastic interaction between holes in 
an elastic sheet was, during the first half of the last cen- 
tury, the subject of much research [9, 14-16, 18, 19]. 
Most of which was dedicated to finding the maximum 
stress felt along the perimeter of each hole. While linear 
elasticity may provide reasonable solutions to such anal- 
ysis, we will demonstrate that linear theory breaks down 
upon further investigation. 

The simplest system accounting for the interaction be- 
tween holes is an infinite elastic sheet containing two 
holes of radius R, whose centers are separated by dis- 
tance 2d. This sheet is then subjected to uniform tension 
P. The analysis of this system will closely follow that of 
Ling [19]. Bipolar coordinates, defined by 



a sinh a 



cosh a — cos P 



y 



a sin/3 



cosh a — cos [3 



(33) 




/3 = 



FIG. 1: Curves of constant a and /3 are circles in bipolar 
coordinates. 



for P £ [0,27r),a G (—00,00), are the natural choice for 
this problem, as lines of constant a or /3 are circles in the 
xy-planc defined by + {y — acot/?)^ = csc^ /3 and 
(x — acotha)^-f — a^csch^a. The system of two equal 
holes corresponds to a = ±s, s = cosh~^(d/i?), and 
a? = cP — R^ . When the system is under uniform tension 
P, the components of the stress tensor are (^xx — i 
Uyy — P, and Uxy — 0. By integrating the Eqs. (21), 
the stress function for an infinite system under uniform 
tension is 



P 



Pa^ cosh( 



cos j3 



2 cosh a — cos (3 



(34) 



Using the method outlined in the previous section, we 
undertake the tedious calculation, detailed in Appendix 
A, to find the equilibrium configuration of this system. 
The results for uniform compression and tension are dis- 
played in FIG. II B 2. Upon further analysis of the com- 
pressed system, large enough values of P yield overlap- 
ping solutions for the displacement vectors, signaling the 
breakdown of the linear theory. 



III. THE LINEAR THEORY OF NONLINEAR 
ELASTICITY 

Even were there no instability in the linear theory of 
two elastic holes, the sheer complexity of the equations 
would make calculations of increasing numbers of holes 
a nearly impossible task, and understanding the mech- 
anism by which the holes collapse and the shapes they 
form requires a nonlinear theory of elasticity. Thus, we 
turn to the theory of cracks for inspiration. In the linear 
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FIG. 2: An clastic sheet with two circular holes cut out (a) 
is subjected to uniform tension (b) and compression (c). The 
dark blue curves arc the the boundary of the holes. The other 
curves show deformations of the circles in (a) to aid the eye. 



theory of elasticity, cracks can be described by a con- 
tinuous distribution of parallel dislocations [8, 20]. The 
stresses in a body due to a crack are the same as the 
stresses in an isotropic body with a distribution of dislo- 
cations with the same height profile as that of the crack. 
As a first approximation, we model each of the collapsed 
holes as a pair of oppositely charged dislocations, known 



as a dislocation dipole [21]. This formalism allows us to 
recover the same physics by the simple numerical mini- 
mization of algebraic equations, once described by a com- 
plex system of coupled differential equations. 



A. Filling a Crack 

Our model system consists of an isotropic solid where 
each collapsed hole is represented by a thin line of mate- 
rial that has been taken out of the system, or a disloca- 
tion dipole. The Burgers vector for a dislocation dipole 





FIG. 3: Two dipoles of strengths di and d2 are separated by 
R. 



of strength b with dipole vector d located at r is 
b (x) = z X db 
or in Fourier space 



{ x-^-rj + (x + ^-r 



(35) 



b{q) = 2i6z X de'«''^°^^sin 



qd 



cos{e - eo) 



ibqzxde*''''^°'^cos{e-eo). 



(36) 



where 9 is angle of q and 9o is the direction of the dipole 
and we have, in the spirit of the dipole approximation, 
taken the lowest order term in d = 

The interaction energy of two dipoles di and d2 both 
of strength b separated by R is given by 



E = 



Y2b^did2 



cPq- 



q X {z X di) ■ q X {z X d2) 



{iq cos{e - ei)) (-ige-'«^^°'^ cos(6» - 62)) , 



(37) 



where di is at the origin. After carrying out the integration (see Appendix B 1), the pairwise interaction between two 
dislocation dipoles is 



Y2 b^rflrfz 



cos(^i + 62) sin 61 sin ^2 + ^ ) • 



(38) 
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Note that the interaction energy is invariant under 9i — > 
Oi + TT and 62—^02+ TT, which reaffirms each collapsed 
hole is represented by a line, not a vector. The total 
interaction energy of an array of dislocation dipoles is 
merely a sum of all pairwise interactions, because we are 
using linear theory. The centers of the initial holes set the 
position of each dislocation dipole, but they are allowed 
to rotate freely. The equilibrium state minimizes the free 
energy over the angle each dipole makes with respect to 
a fixed axis. 

I 



over each of the angles, which arc measured with respect 
to the X-axis. Minimizing with respect to the four angles 
we find 61 = 04 and 62 = O3 by symmetry and 

sin 46'i - cos 26*1 -4 sin 2 (6*1 +02 ) = 

sin402 + cos26i2 -4sin2(6ii +6'2) = 0. (40) 

The minimum is 

While one might have postulated that the lowest energy 
configuration would have = and 62 = 7r/2, it turns 
out that the energy is slightly lowered if these angles 
are slightly shifted. This is due to the finite size of our 
system - as we shall see, for larger systems the dipoles 
align along the crystal axes and there is a boundary effect 
which distorts the dipole directions at the edges. 

2. We Had to Count Them All: n x n Systems of Holes 

One might wonder how we can study larger and larger 

systems since the interaction only falls off as Be- 
cause the interactions is between dipoles, the interaction 
energy at large distances decreases becaiisc the angle of 
the dipoles rotates around the circle. Indeed, consider the 
interaction of a single dipole at the origin with N"^ — 1 
other dipoles in an Na x Na lattice. The angular depen- 
dence in (38) will wash out the power law if the dipoles 
rotate through 27r uniformly. As we will see in the next 
section, this is precisely what happens as shown in Fig. 
5. 

For these larger arrays of holes, the sheer number of 

coupled equations makes it impractical to find solutions 
by hand, and we turn to numerical methods to find the 



1. The 2x2 Diamond Plate Plaquette 



The simplest case consists of four holes located at 
{±a/2,±a/2}. Because all four holes have the same ra- 
dius, their dipole vectors should have the same magni- 
tude. We minimize the energy functional, composed of 
the sum of six pairwise terms. 



f)sin(03-f) 

(39) 



orientations of the ground state of each lattice. We find 
that the diamond plate order of the 2x2 plaquette per- 
sists for larger and larger samples with increasing align- 
ment along the lattice directions. 

I - 
- I 

FIG. 4: The groundstate orientation of the 2x2 plaquette. 

B. Stretching the Lattice 

In experimental systems [4], the elastic sheet was 
stretched in a specific direction before allowing it to 
swell, leading to a background stress axx = Tcos(j), 
(Tyy = T sin (f). To calculate the coupling energy between 
the stretching and the dipole angle, we will rotate the 
system, such that the x-axis is defined by the direction 
of stretching, or Gxx = T and the dipole is located at the 
origin, or 

b{q) = ibqzxdcos{e- Bo + (!)). (42) 

It is most appropriate to use the energy functional E = 
2^ / (Px (V^x) . In Fourier space, the dipole term is 
given by V^x = i^ik^bi{q) [21] and the stretching term 



-£'2x2 = — ^ rcos(0i + ^2) sin ^1 sin 62 + cos(03 + 9i) sin ^3 sin 64 + cos(0i + ^3 — tt) sin(0i — 

+ cos(6i2 + 6I4 - tt) sin(6l2 - f ) sin(6'4 - f ) + | cos(6'i + 6*4 - f ) sin(6ii - f ) sin(04 - f ) 
+ i cos(02 + ^3 - f ) sin(e2 - ^) sin(03 - ^)] , 

I 
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(c)<^lft^Ib±ciTV^^^(d)^a^|^^ 




FIG. 5: The groundstate orientations for square lattices of 
10 X 10 (a), 13 X 13 (b), 17 x 17 (c), and 20 x 20 (d). 
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by V^x = T^S{e). The energy is 



1 
2Y2 



-T 



Q 



bdcos^ie-Oo 



S{6) +bdcos^{e -6Q + (f>) 



(43) 



from which the couphng energy is 



-^stretch — 



Tbd 
Tbd 



J qdq J d9^S{e)cos'^{e-eQ + (f)) 
cos^(j>-eo) (44) 



The new term causing the dislocation dipolcs to align 
with the direction of stretching competes with original 
interaction energy, favoring the diamond plate pattern. 
Following the same minimization procedure as before, we 
find that for small tensions the diamond plate pattern is 
only slightly perturbed, and for large tensions, the dis- 
location dipoles align along the direction of stretching, 
which may be seen in Fig. 6. 



IV. CONCLUSION 

We have created a model system for an elastic sheet 
with a square lattice circular holes cut out of it. When 
the sheet is swollen, or, equivalently, subjected to uniform 
tension, the holes snap shut. Their major axes align into 
a diamond plate pattern with long ranged order. While 
other methods of calculating the orientational order of 
the holes rely upon nonlinear elasticity theory and finite 



FIG. 6: The 4x4 and 7x7 system of holes are stretched 
by tension Tx. For small tensions (in (a) and (d) T = 
0.1^^^^), the diamond plate groundstate is only slightly 
perturbed. Whereas for large tensions (in (b), (c), (e), and 
(f) T = 2^^-^), the holes align along the i-axis, the direc- 
tion of stretching. In (b) and (e) the angle of stretching is 
^Stretch = 7r/4 from horizontal. 



element simulations, we use simple linear elasticity the- 
ory to obtain the same results. Our system also corrob- 
orates experimental results of the sheet under external 
forces. This method may easily extended to holes on 
other lattices. It is difficult to extend it to an infinite 
lattice, since the minimization would then be over an in- 
finite number of angles. An Ewald type summation may 
be used for an infinite system whose unit cell is the 2x2 
plaquctte. However, this is unlikely to lead to new in- 
sight, as the diamond plate order is clearly maintained 
for large systems. Moreover, as we show in the appendix, 
the dipole interactions are the dominant terms even for 
more general elliptical holes. 
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APPENDIX A: LINEAR ELASTICITY IN 
BIPOLAR COORDINATES 



from (TaQ+(T/33 and (t,,.,, — (Tari, respectively. There exists a 
function G which satisfies equation (A6) for which ^ = 

^ - ^huo, and ig = -0 + ^hup. Using these 
equations, we eliminate and Uf) from equation (A5) 



1. The Basics 

For clarity's sake, we take a moment to explicitly write 
out the useful equations in bipolar coordinates. From 
equations (8), we simply read off the components of the 
strain tensor 



1 dUa 

h da 



1 dh 



1 du/3 1 dh 



2u. 



a(3 



da\h)dp\h 



h? da 
(Al) 



The stress tensor may be written as a function of the 
Airy stress function x, from equation (23) 



ld_ 
hd/3 
ld_ 

h da 



Idx 
hdf3 

Idx 
h da 

,2 



+ 



1 / d^x 
h? \dadl3 
d_ 

da 



1 dh dx 
h^ da da ' 

1 dh dx 

h^'m'dd 

1 dh dx 1 dh dx 
hd/3 da h da d(3 

d f 1 dx 
dp [l^d^ 



(A2) 



1^' 
h^ dp, 

or in terms of the relative displacement vector, from Eq. 
(17) 



0-/3/3 



Y2 



1-C72 

Y2 

I- a 



{Uaa + CTW/3/3) , 
2 {Upi3 + aUact) , 



1 + a 



Ua(j. 



(A3) 



At first glance, it seems an insurmountable goal to solve 
the differential equations for the displacements in terms 

of the Airy stress function. However, following [14], they 
become much more tractable if one considers the terms 



Yo 1 



1 — (J h^ \ da 

_ 1 / d'^x (^"^X 
h? \d^ ^WJ' 



d{hua) _|_ d{hufj)^ 



dp J 



(A4) 



Y2 



1 



d /Ua 

da \ h 



d nip_ 
dp\h 



da \h^da) dp WdpJ ^ ' 



or, equivalently, we arrange equations to obtain 

9_(dx_ Y2{hua) \ d(dx_ ^2(^^/3) 
da\da 1-a J dp \dp 1-a 



-0, 
(A6) 



d^ fG 
dadp \ h 



{aca -(^00), (A7) 



the left hand side of this equation may be written in this 
manner because g^g^^ = 0. Thus, the relative displace- 
ment vectors are given by 



U0 



I ■ /dx _ dG 
Y2 h \d^ 'dp 
1-a l /dx dG 

Y2 h \dp ^ 



(A8) 



2. Fixing a Hole (Outtake) 



Recall, X is the solution to the biharmonic equation 
V^X = with respect to the symmetries in our system. 
The biharmonic equation when written in terms of the 
function x/h has the simplified form, 



94 



94 



(— — 

\da^ ^ dp^ ^ da^dp^ 



92 



d^ 



^da^ dp^ 



' (A9) 

Similarly, the components of the stress tensor are given 

by, 



fa/3 



1 d^ 
hdp^ ~ 
1 d^ 
h da'^ 
1 92 



hdadp 



sinh a d sin P d cosh a 

a da 
sinh a d 

a da 

©■ 



a dp a 
sin/3 d cos/3 
a dp a 



X 
h 

(AlO) 



Our system, while undergoing uniform hydrostatic 
compression, is described by the stress function in equa- 
tion (34), 



y(0) 

— — = -— (cosh a + cos P) , 

n Zi 

or, equivalently, is given by the components of the stress 
tensor 



.(0) 



JO) 

^00 



(All) 



We aim to find solutions to the biharmonic equation that 
are even in both a and p. Thus, the Airy stress function 
is given by 
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/(I) 



= C (cosh a — cos (3) log (cosh a — cos /^) + / ^ 0n (o^) cos ^/^; 



(A12) 



n=l 



where 4>n{c() = An cosh(ri + l)a + B„ cosh(n — l)a. The components of the stress tensor corresponding to this Airy 
stress function are 



r(l) 



C 1 

(cosh2Q! — 2 cosh a cos /3 + cos 2/3) H — (/'ifa) 

2a a 



1 °° 1 

+ — V - [/„+i(a) - 2cosha/„(Q!) + /„-i(a) - 2 sinh a5„ (a)] cosn/3 
2a ^ — ' " 

C 



2a ^ n 

n— 1 

(cosh 2a — 2 cosh a cos /3 + cos 2/3) + —(j)i{a) — t— (/'"(a) 



(A13) 



- ^ [<+i(a!) -2 cosh a(f)n{a) + <_i(q!) + (n + 2)0„+i(a) + 2sinha0;(a) + {n- 2)0„_i(a)] cosn/3 (A14) 

n=l 

CT^^j = sinhasin/3 - — ^ [fif„+i(a) - 2coshafif„(a) +5„_i(Q!)]sinn/3, (A15) 



where /n(a) = (n+l)n(n— l)(/i„(a) and S(„(q!) = n0^(a). 
The boundary conditions require there be no stress at 
infinity (a = 0), hence 



J2 [An + Bn] = 0, 



(A16) 



n=l 



and the normal and tangential stresses must vanish along 
the edges of the holes, located at a = ±s. Thus, the 
constants An, _B„ and C must satisfy the following re- 
currence relations for n > 2 



fn+iis) - 2 cosh s/„(s) + /„-i(s) 
gn+i{s) - 2coshs9„(s) +gn-i{s) 



2sinhs5„(s)j(A17) 
0, (A18) 



subject to the conditions 

/2(s) - 2sinhs5fi(s) 
s-2coshs/2(s) + /3(s) 
2coshsfifi(s) -52(5) 



-2- -C cosh 2s, (A19) 

2 cosh s, (A20) 

2sinhs5f2(s) + 2C, (A21) 

2Csinhs. (A22) 



Using equation (A18), we find gnis) = ciA" + C2A2, 
where Ai and A2 are roots of the characteristic poly- 



nomial - 2cosh.st" + f'-^ = 0, yielding ,9„(,s) 



cie 



+ C2e" 



However, the stress must be finite 



everywhere, thus, C2 = 0. Equation (A22) com- 
pletes the relation for .g„(s) as ci = 2Csinhs. Be- 
cause, equation (A17) is a non-linear recurrence rela- 
tion, we must consider fn+2 — 2coshs/„+is -|- /„s — 
^(/n+i(s)-2coshs/„(5) + /„_i(s)) = 0. Two of 
the roots of the characteristic polynomial for this equa- 
tion are degenerate, the equation is /«(«) = rfie""* -|- 



done 



d^e^" . From the boundary conditions, wc de- 



termine fn{s) = — 2ii' (coshs -|- nsinhs) e""*^. From the 
definitions of fn{(^) and 5„(a), we find that the coeffi- 
cients An, Bn and C satsify. 



An 

Br, 



2C 



' sinh J 



' sinh ; 



n{n + 1) (sinh 2ns -|- n sinh 2s) ' 



-2C 



e sinh ns + ne^ sinh s 
n{n — 1) (sinh 2ns -|- n sinh 2s) ' 
C 



with Bi = — tanh s cosh 2s -|- P, 



(A23) 



and 



C = -pi - -htanhssinh2s-4^ 



n=2 



e sinh ns -|- n sinh s (n sinh s -|- cosh s) 



n{'n? — 1) (sinh 2ns -|- n sinh 2s) 



(A24) 



Now that we have equations for the stresses everywhere, we may now solve for the field of relative displacement 
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vectors. Recall, our function G is given by, 



2h 
1 + <J 



2C i tan ^ Ttanh ^ cot ^ \ cos 



3 /? + tan ^ f coth ^ tan ^ j cosh a j + (oi) sin n/3 

^ ' ' n=l 



(A25) 



where X = X^"^ + X*-^-* and '(/'(o^) = sinh(n + l)a + -B„sinh(n — l)a. To calculate the displacement field, we will 
need to know the following relations. 



1^ = /i (sinha(-P/icos/3 + C) - ^ 



9x 

dG 

da 

dG 
dp 



n=l 

oo 



—(f>n{a) sinha — (/"^(a) 



cosn/3 



= h [ sin p{—Ph cosh a + C) — ^„ (a) 



1 + a 

2/i 
'l + a 



n=l 



■ sin /3 cos n(5 + n sin n/3 



"^h l^f h . ^ ^ . \ ^ 
C I 77— smh a cos p + sin p I + > 

^ ^ n—l 



-tpn{oi) sinha — i^'nia) 



h 



C ( 77— cosh a sin /? — sinha J + ^ 



sin n/3 



■ sin /? sin n(3 — n cos n/3 



(A26) 
(A27) 
(A28) 
(A29) 



r 



These equations, together with Eqs. (A8), complete our 
description of the system of two holes under hydrostatic 
compression, which may be seen in Fig. 2. It should 
be noted that in the compressed system there is an in- 
stability for large enough values of P wherein the dis- 
placements intersect each other, causing overlap in the 
system. 

Due to the nature of problems in the theory of lin- 
ear elasticity theory, the most useful identities involve 
the directional cosines relating the {a,/3} to the {x,y} 
coordinates. The unit vectors in the new system are 
given by d = ^ (If A + l^y) and ^ = ^ (§& + %y) . 

The orthogonality condition a. ■ ^ = Q implies + 

0. The directional cosines are related because 



9y dy 
da 83 



K 



(if) (l+(£/|f) ) = (t) (! + (-§/§) ) 



1 dx 
ha da 



J_dy_ 

h(}dp' 



J_dy_ 

ha da 



J_dx 
'h~pW 



(A30) 



where the sign is chosen such that both {x,y,z} and 



{q;,/3, i} form right-handed orthonormal triads. Thus, 
under the change of coordinates {x,y,z} {a,P,z}, 
a rank-2 tensor Aij — UijXiXj transforms as A^i, = 

o-ixi^^i^^i^ = aij {^^ ■ ^i) (^L ■ ^j) i^-, which may be 
written as 



^111/ — 



1 dxi dxj 
h^K d^i^ d^^ 



Ctii • 



(A31) 



The derivatives of the directional cosines can be made 
from linear combinations of derivatives of the orthogonal- 
ity condition -i^^ + -f^-l^ = and the scale functions. 

•' da dp ^"P^ a X 

For example, to find ^(7^5^); first note there are two 
ways of obtaining this derivative; directly. 



d_ 

da 



1 dx 
ha da 



1 dha dx 1 d'^x 

1 . 

da da ha do? 



(A32) 



and by taking the derivative of the product of a scale 
function and orthogonality condition with respect to a 



d ( \ dx\dx ^ f ^ ^^<^ ^^y ^y ^^-^ ^^y ^y\ 

da \hadaj db ha \ ha da da dp do? dp dadp da dadp da J 

^ J_dh^dy_dy_ _ l_d^dy_ _dha . 
hi da da dp ha da^ dp dp ' ^ ' 

where we have made use of the definition = ^ + W~d8 • ^^^t multiplying Eqn. A32 by fe and Eqn. A33 
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by (If )^ and taking their sum, this becomes 



1 dx\ 1 dha dy / dy dx dy dx \ dha dx 1 dy ( dy d^x dx d^y \ 

1 dha dy hfj dha dx 1 dy hfj dha 
da d(3 ha d(3 dfi ha dj3 ha da 

By following a the same procedure, the formulae for the derivatives of the directional cosines are 

d f 1 dxi \ 1 dha dxi d f 1 dxi \ 1 dhf) dxi 



(A34) 



da \ha da J h"^ d(3 d(3 ' dfi \ha da J hah/} da dj3 ' 

d f 1 dxi \ 1 dha dxi 9/1 dxi \ 1 dhfj dxi 



da \hi3 d(3 J hahp d(3 da ' d(3 \h/3 dp J h^ da da 



(A35) 



where « = 1,2 and xi = x and X2 = y- Note that 
the cross partial derivatives of the directional cosines are 
equal. 



d / d / 1 dxi 
dp \da \ha da 

d / d / 1 dxi 
'dp \d^ \Mi~dp 



d_ 

da 
d_ 

da 



d_ 

'dp 
d_ 
'dp 



1 dxi 

ha da 
1 dx.^ 
hp dp 



which leads to our final identity 



d_ 
'dp 



1 dha 
hp dp 



d_ 

da 



]_dhp 
ha da 



0. 



= 0, 

= 0(A36) 
(A37) 



3. Identities in Orthogonal Coordinates 

We include these for completeness. These were neces- 
sary for us to study the Airy stress formalism in orthog- 
onal coordinates and we did not find these, presumably 
known identities, in any reference. 

Consider the general set of orthogonal coordinates, 
{a{x,y),P{x,y)}. The new basis preserves length of the 
differential line element, ds^ = dx^ + dy^ = h^^da^ + 
= + %dvf + hlifjx + §dyf, defin- 

ing the scale functions = (|^) + (^) and h^^ = 

+ (f^)^- However, it is often more useful to con- 
sider Cartesian coordinates as functions of the new or- 
thogonal ones, {x{a, P),y{a, P)}, which yield an equiva- 



lent statement of the scale functions = (ff ) + (ff 
and hp = (f|)^ + (|^)^-By transforming from or- 
thogonal back to Cartesian coordinates, the differen- 
tial line element gives + 

set of identities: 



1 and 

1; from which, we obtain the first 



dx ^2 (^'^ ^2 
da " ' da " dy 



(A38) 



dp ^dx' dp ""dy' 

APPENDIX B: FILLING THE CRACKS: THE 
DISLOCATION DIPOLE 



We devote this Appendix to the technical details of 
the mathematical manipulation required to compute the 
energetics of a lattice of dislocation dipoles. 



1. The Dipole Term 

Despite the complex form of the integral in the dislo- 
cation dipole interaction energy, when completed it has a 
surprisingly simple form. The integral becomes tractable 
by transforming to polar coordinates, and then manipu- 
lating the trigonometric functions. In polar coordinates, 
the integral in Eq. (37) becomes. 



^ Y2b''dld2 p . . 

Employing Bessel function and trigonometric identities, 
we find 

E = -^^^^ (^cos(^i+02)sin0isin02 + ^) ■ 

(B2) 



qR cos(6»)) cos2(6» - 6»i) cos2(6» - 62) (Bl) 

I 

2. Why We Can Ignore Higher Order Terms 

Our goal is to prove that we need only consider the 
first order dipole-dipole term when considering collapsed 
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ri(.T) 



(b{x) 



d/2 



holes. For simplicity sake, we will study an elongated 
shape that is symmetric about both its semi-major and 
lix)// semi-minor axes, (see FIG. 7). Here, we consider only 
shapes for which the ratio of minor and major axes 
2a I d <C 1. In the theory of cracks, the height profile 
of a crack h{x) may be constructed from a continuous 
distribution of finite parallel edge dislocations. A finite 
edge dislocation of length I can be thought of as a two 
infinite edge dislocations of the same strength but oppo- 
site charge which is given by Burgers vectors situated at 
6+ = b5{x—t/2) and6_ = — 6 ^(a;+.^/2), in other words, 
a dislocation dipole. Given the height profile of a shape, 
it is trivial to construct it from such dislocation dipoles, 
h{x) = d{x). Thus, the interaction energy between such 
a shape made from dislocation dipoles and a single dis- 
location dipole d2 located a distance Rq away is given 
by: 



-d/2 



FIG. 7: Set up for dipole expansion of an elongated shape 
symmetric about both the x- and y-ax.es, in this case, an 
ellipse. 



E = -^^V^ /" (cos (cl^ix) + 02ix)) sin {cPix)) sin {O^ix)) + . (B3) 

By repeated application of the law of cosines, the functional forms of R{x), (l){x) = 7r/2 — a{x), and 02{x) are: 

R'^{x) = Rl+x'^ -2 Rqx cos (^^+eoJ = Rl + x^ +2 Rox sin 6*0 (B4) 

, . -if x + RQsm9Q \ 
a(x) = cos — (B5) 

\^Rl+x^ + 2Roxsm9o J 

02{x) = e2,-0o-a{x)+n/2. (B6) 

The energy density is 

^ _y2^_d2 (_ cosW - 2 a(x)) cos(a(x)) cos(^ - a(x)) + \] , (B7) 

TT i?^-|-a;2 + 2i?oa;sin^o V \ J/ \ \ J/ v ^ " 4 J ' ^ ' 



where 6 = 62^ — 6q- This may be vastly simplified by expanding the angular terms: 

cos6{x + Ro sin OofR^ cos^ 9^ + 2 sin 9{x + Rq sin ^0) Ra cos ^0 
Rl+x"^ + 2Rqx sin 6*0 



~ K^yjawyuj -r J.>A}aii.i.wQ J i,q v^vjo wij -r ^ 01.1.1. \j yj, -r J. >A} alii, wij j K^yja wij /no\ 

co&{9-2a{x)) = d2 , ^2 , o D (B8) 



,7: , cose'(x + i?osin0o) +sin^i?ocos6'o 

cos{9-a{x)) = . , • (B9) 

\/Rq + x'^ + 2Rqx sm t'o 
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The energy density becomes: 

^ Y2bU2 d{x) nRl+x^+2RoxsmGof ^ (x + A)A 

= 3 y c„ y (BIO) 

where A = i?osin6'o and B = RqcosOq, and the coefScients c„ are given by cq = Ci = -B^ cos^ sin^, C2 = 

B^{cos^6 — 2 sin^ (?) + C3 = — 3i? cos9 sinO, and C4 = — cos^ 6 + 1/4. While the above energy is for general 

shape of dislocations, we choose a shape to do the actual calc ulation. For simplicity sake, we choose an ellipse of 
major axis do and minor axis 2a. Thus, d{x) = do-\/l — x'^ jc?. With the change of variables, y = x/a, our energy 
integral becomes: 

E=-^^^ t dy WlIZ sEc^E „ A"-(a,)-. (Bll) 

-^0^ i-i '(l + (^,)2 + 24,sin^o) "^o^i'o-Kn-m)! 

Expanding the denominator for ^ 1, we find, ^l + (-^y)^ + 2-^2/sin(?o) ~ 

(^1 - 6-|^sin(?oy + 3(-|^)^(8 sin^e'o - 1)?/^ H )• Now we need only do the integral /j^^ dyi/T^-y^y" = 

2(1 — (— 1)") /o^ y'l — y'^y^ , for, since \/l — y"^ is even, this integral is zero for odd integer n This is very 
simple using beta functions, which we may see by the change of variables t = y'^: 

2 t dy'^/l^y"" = /'^x/r^i"/2= rdft("+i)/2-i(l-f)3/2-i 
Jo Jo vt Jo 

_ p/l + » 3y r(l±i^)r(|) _ 7r(l + (-l)»)(n-l)!! 

- V 2 '2^ r(2+f) 2t+i (f + 1)! ^ ' 

The interaction energy is thus: 

^ n,=n m.=n ^ ^ \ \ Z 



While this appears to be a complicated expression, let us, for the moment consider only the first order term in a: 

Y2b''d2doaJ-. ^„ Y2b'^d2doa /B^ ^ . ~ ~ 

U n=0 ^ ^ 

+ A^B^i cos2 6i-2sin^^ + ^ ) - 3 Bcos6'sin6» + A* ( t - cos"" 



Y2 6^ d2doa f cos'' ^0,/ sn ■ a -3 ■ 3 n a\ 2 ■ a 

— o h I cos 6^0 sm fc'o — 3 sm Oq cos ) cos sin 6/ 

4i?n V 4 ^ ^ 



1-0 

+ cos^ do sin^ (^cos^ 0-2 sin^ 61+^^+ sin" 9o (^^ - cos^ 0j j . (B14) 

This does not appear to have the same functional form as the original energy. However, with the help of some 
trigonometric identities, we begin to simplify the energy: 

cos" 00 + sin" 6I0 (1 - 4 cos^ ^) = \ ^in^ 26^0 cos 2^ + 1 - 4 sin^ 9o cos^ (B15) 

( cos^ sin Oq - 3 sin^ Oq cos ^0) cos sin 61 = ^ (2 cos 26*0-1) sin 26*0 sin 20 (B16) 

(cos^6i- 2sin2 6i+ cos^ 6*0 sin^ 6*0 = ^ cos 26isin2 2610. (B17) 
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We now use these identities in the energy to find: 

I2 C^O 0^ / 1 . 2 o/i r,;; 1 -2/1 2r\ ^ 



= .r)2 \ o sin 26*0 cos 26* + - - sin^ Bq cos ^ + - (2 cos 2^0 — 1) sin 2^o sin 2^ + - cos 26'sin^ 2Bq 

4 /in 



Y2b'^d2doa fl . 



- sin 2^0 (sin 29o cos 26^ + cos 26o sin 261^ - sin^ ^0 cos^ ^ - sin 260 sin 2^ + i j . (B18) 



Now we substitute = — Oq to find, 



Ynb do dna / 1 „ „ o„ o ^ „ „ 7: ~ 1 



Eq = —2 - sin 29o sin 2^2o ~ sin ^0 cos 6 — sin ^0 cos ^0 sin 9 cos + 7 

4 itg V 2 4 



F2 6^ d2 do a f^^^^^^ ^ ^^^^ ^.^ ^.^ 02, + ]]. (B19) 



4i?2 y—^-- ■ "-0. ^0 . 4 



I 

So we have found that to first order in a/i?o an elhpse of Higher order terms contain more powers of a/Ro, which 

major axis d and minor axis 2a interacts with a disloca- can be neglected for very thin shapes or for dipoles that 

tion dipolc of strength ^2 a distance Rq away like a dislo- are very far away, 
cation dipole with effective dipole strength d = -k doa/A. 
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